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Abstract. An algebraic multilevel iteration method for solving system of linear algebraic equations arising 
in H{Q, curl) space is presented. The algorithm is developed for the discrete problem obtained by using the 
space of lowest order edge elements. The theoretical analysis of the method is based only on some algebraic 
sequences and generalized eigenvalues of local (element-wise) problems. Explicit recursion formulae are 
derived to compute the element matrices and the constant 7 (which measures the quality of the space 
splitting) at any given level. It is proved that the proposed method is robust with respect to the problem 
parameters, and is of optimal order complexity. Supporting numerical results, including the case when the 
parameters have jumps, are also presented. 



1. Introduction 

Consider the finite element discretization of variational problems related to the bihnear fonii 
00; (1) 9<{u,v) := a{u,v) + ^{cmlu,cmlv), a,ySeR+, 

defined on the two-dimensional Hilbert space 
<[^'; (2) H{a,cml) := {v e {L^{Q.))^ : curl v e L^{Q.)}. 

^ I Here Q is a Lipschitz domain, curl v = dxV2 — dyV\ is the scalar curl of a two-dimensional vector 

' V = \v\,V2Y 1 and (•, •) denotes the inner-product in L^(Q). For a = p = 1, the bilinear form ([T]i is 

"j^ ■ precisely the inner-product in H{Q., curl ). 

Associated with the inner-product 1-1, there exists a linear operator A : = al + yScurl curl, which is 
determined by the relation 

(M ; (3) {Au,v) ='H{u,v), ^veH{a,cm\), 

\ where the vector curl of a scalar function w is defined as curl w = \dyW, —dxwY- Given a finite element 

\C> ' space "Vh of //(Q, curl), the symmetric and positive-definite (SPD) operator A/, : 'Vh "Vh, which is 

^ . the discretization of the operator A together with natural boundary conditions, is determined by 

(4) {AhUh,Vh) ='H{uh,Vh), Vv/, eO//,. 

^ \ The operator equation Au = f, for / e (L^(Q))^, then leads to the following discrete problem 

(5) Ai,Uh = fh, 

which is uniquely solvable. Such problems frequently occur in various contexts in electromagnetism, 
^ I e.g., low-frequency time-harmonic Maxwell equations E4\ . or some formulations of the (Navier -) 

^ ' Stokes equations |[T31l . Therefore, developing fast solvers for large system of equations ([5) is of sig- 

- ■ ■ nificant importance. 

Preconditioning methods for such linear systems within the framework of domain decomposition 
(overlapping Schwarz), multigrid and auxiliary space methods have been proposed by several authors. 
The first results for multigrid in H{D., curl ) (within the framework of overlapping Schwarz methods) were 
obtained by Hiptmair in |[T4l . A unified treatment of multigrid methods for //(H, curl) and //(Q, div) 
was presented by Hiptmair and Toselli in 1151 . However, the condition number estimates of their precon- 
ditioned system were not robust with respect to the parameters a and j3. Arnold et al. ||T] employed the 
multigrid framework by developing necessary estimates for mixed finite element methods (FEM) based 
on discretizations of H{Q.,div) and //(Q, curl), and thereby obtained parameter independent condition 
number estimates of the preconditioned system. Pasciak and Zhao studied the overlapping Schwarz 
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methods for //(Q, curl) in polyhedral domains in E51 . and Reitzinger and Schoeberl studied algebraic 
multigrid methods for edge elements in Auxiliary space preconditioning, proposed by Xu in 
was studied for //o(^>curl ) (the space //(O, curl) with zero tangential trace) by Hiptmair et al. [i6i|. 
Nodal auxiliary space preconditioning in //(O, curl) and //(Q,div) was studied by Hiptmair and Xu in 
ifTTl . and the proposed preconditioner was robust with respect to the parameters a and /?. 

In this paper we develop optimal order iterative solver for the solution of the linear system ([5) based on 
algebraic multilevel iterative (AMLI) methods. AMLI methods were introduced by Axelsson and Vas- 
silevski in a series of papers ||3]|4l|5]|6]- The AMLI methods, which are recursive extensions of two-level 
multigrid methods for FEM f2\, have been extensively analyzed in the context of conforming and non- 
conforming FEM (including discontinuous Galerkin methods), see |l8l[T2l[T8j|20l|2TJ|22l. For a detailed 
systematic exposition of AMLI methods, see the monographs IIT91I311 . To reduce the overall complexity 
of AMLI methods (to achieve optimal computational complexity), various stabilization techniques can be 
used. In the original work EH], the stabilization was achieved by employing properly shifted and scaled 
Chebyshev polynomials. This approach requires the computation of polynomial coefficients which de- 
pends on the bounds of the eigenvalues of the preconditioned system. Alternatively, some inner iterations 
at certain (or all) coarse levels can be used to stabilize the outer iterations, which lead to parameter-free 
AMLI methods 15] |6l [181 125] |26l [271. These methods utilize a sequence of coarse-grid problems that 
aie obtained from repeated application of a natural (and simple) hierarchical basis transformation, which 
is computationally advantageous. Moreover, the underlying technique of these methods often requires 
only a few minor adjustments (mainly two-level hierarchical basis transformation) even if the underlying 
problem changes significantly. 

Our analysis is based only on some algebraic sequences and the generalized eigenvalues of local 
(element-wise) problems. We derive explicit recursion formulae to compute the element matrices and 
the constant y (which measures the quality of the space splitting) at any given level. The method is shown 
to be robust with respect to the parameters, i.e., the results hold unifomily for < a,/3 < az>. 

The remainder of this paper is organized as follows. In Section |2] we briefly discuss the finite element 
discretization of the model problem ([T|l using the lowest-order Nedelec space. Section [3] starts with a 
brief description of the AMLI procedure (in Section 13. It . and is followed by the construction of the 
hierarchical splitting of the lowest-order Nedelec space (in Section !?^ . In Section [X?] a local two- and 
multi-level analysis is then presented and the main result is proved. Finally, in Section H] we present 
numerical experiments. These include the cases with known analytical solution (a = yS = 1), fixing one 
of the parameters a or jS and varying other from 10^^ to 10^, and jumping coefficients. The conclusions 
are drawn in Section |5] 



2. Finite element discretization using Nedelec elements 

We consider the tessellation of Q using rectangular elements, and choose the reference element 
^ as [—1,1] X [—1,1]. Let P^i ,r2 {K) denote the space of polynomials of degree ^ ri in ;c and ^ r2 in y. 
Also, let Pr{dK) denote the space of polynomials of degree ^ r on dK. For the construction of "V/,, we 
use the space of lowest-order edge elements (Nedelec space of first kind), which is denoted by N°. The 
space N°(^) is defined as 

(6) ^\k) = PoAK) X Pi^K) = f^v{x,y) = 

Thus, the local basis for has dimension 4. Moreover, for vq e N"(^) we have 

(7) curi^;oe^o,o, vq ■ t\g^ e Po{dK), 

where t denotes the unit tangential vector to the element boundaries. For further details the reader is 
referred to, e.g., |[24l . 

Now let F : ^ ^ be a diffeomorphism of the reference element K onto a physical element K, i.e., 
K = F{K). By J' we denote the Jacobian matrix of the mapping, and by J'o its determinant, which are 
defined as 

1^)' ^D = \detj\=d,xdyy-dyxd,y>0. 
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Then we have the following transformation relations: 



(8) 



w = J' w; curl w = J'j^ ^cuii w, \fw e H{K, curl), w e H{K, curl). 



The vector transformation w J'^ w is called the co variant transformation, and w J'^ J'w is the 
well known Piola transformation. 

We denote the element matrix for S^^^u ■ vhy Lk, and for curl u curl v by Ck- For the N" space 
based on uniform mesh composed of square elements, the element matrices Lk and Ck have the following 
structure 



(9) 



1 



2 
1 










1 



The overall element matrix Ak, which is defined as Ak 

2ah^ + 6/3 ah^ - 6/3 

1 ryh^ - 

(10) Ak 



oCk 



1 

6^ 



PBk, is thus given by 
6/3 



-6/3 
6/3 



2ah^ + 6/3 
6/3 
-6p 



^6/3 
6/3 -6/3 

lah^ + 6/3 ah^ - 6/3 

ah^ - 6/3 2ah^ + 6/3 



Letting e = kH^, with k = a/p, the element matrix can be written as 



(11) 



2e 

e - 



e - 
2e 



2e 
e - 



e - 
2e 



Clearly, for all a,p e ]R+, and thus k e we have e > 0. Note that for fixed k, and h —>^ 0, the element 
matrix Ak is dominated by the matrix Ck (which has a non-zero kernel), whereas for moderate values of 
h it is a regular matrix. The near-nullspace of the matrix Ak is given by the nullspace of the matrix Ck, 
which is associated with the local bilinear form Ck{u,v) := (curl it, curl v)k- As we shall see in the 
analysis, the proposed method is of optimal order for all < a,p < oo. 

Lemma 2.1. (Near-nullspace of matrix Ak)- The element matrix Ak (flOl ) is symmetric positive definite 
(SPD). Moreover, the nullspace of the matrix CKfor a general element K with nodal coordinates (x,-, j, ), 
/ e { 1, 2, 3, 4} is given by 

(12) ker(Cif) = span{(l, 1,0, 0)^,(0, 0,1, 1)^,(^1, X2,B,3'4)^}. 

Furthermore, in case of a uniform mesh composed of square elements, the matrix Ck is same for each 
element K and its nullspace is given by 

ker(Cif) = span{(l, 1,0,0)^, (0,0, 1, 1)^, (-1,0,0, 1)^}. 

Proof. Since the coefficients a and /? in (ITOl) are positive, it follows from equation ([5) that Ak is SPD for 
a general element K. Moreover, it can easily be seen that there exist positive constants c\ and C2 such 
that the following inequalities hold: 

(13) ci ^(w, lu) w^Ca:w ^ C2 ^(w. If). 

Here is a function that is linear in one component and constant in the other component, and the 
function ^(s •) is defined by 

— dyW[)^dK. 



^{w,w) := (5vi 
Jk 



Thus (for the paiticulai^ choice of w), we get C^rw = if and only if id = (ai + a^y, 02 + ajxY , for all 
constants o-,, / = 1, 2, 3. By setting one of these constants to 1 and the other constants to 0, we obtain the 
three linear independent vectors {x\,X2,y3,y4Y , (1, 1,0, 0)^, and (0, 0, 1, 1)^, that span the kernel of Ck- 
The remaining part related to unifom mesh is trivial. Note that, in case of a unifom mesh composed of 
square N" elements, since the vector (1,— 1,-1, 1)^ is orthogonal to the kernel of Ck, it is cleai^ that the 
rank-one matrix Ck is of the form c ■ (1, —1,-1, 1)^ ■ (1, —1,-1, 1), for some constant c. □ 
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Remark 2.2. When using the lowest order Nedelec elements, the matrix Ck is always of rank one. In the 
global assembly this yields a matrix C whose rank equals the number of elements in the mesh. That is, 
the kernel of the global matrix C has dimension dim(ker(C)) = he — riK, where he denotes the number 
of faces and hk the number of elements in the finite element mesh. Thereby, the dimension of the kernel 
is slightly more than half of the total number of degrees of freedom. 



3. Algebraic multilevel iteration 

For the solution of the linear system arising from ©, we describe and analyze the AMLI method in 
the remainder of this section. Our presentation follows Reference Il22l . 

3. 1. The AMLI procedure. In what follows we will denote by M'^'^^ a preconditioner for a finite element 
(stiffness) matrix A^^) corresponding to a ^ times refined mesh (0 ^ ^ ^ L). We will also make use of 
the corresponding level hierarchical matrix A^^\ which is related to A^^^ via a two-level hierarchical 
basis (HB) transformation i.e., 

(14) a(^) =7Wa(^)(/W)^. 

The transformation matrix j'^^'^ specifies the space splitting, which will be described in detail in Sec- 
tion |3]2] By A^l^ and A^l'^ , 1 ^ /, 7 ^ 2, we denote the blocks of A^^^ and A^^) that correspond to the 
fine-coarse partitioning of degrees of freedom (DOF) where the DOF associated with the coarse mesh 
ai^e numbered last. 

The aim is to build a multilevel preconditioner M^^^ for the coefficient matrix A^^^ := A/, at the level 
of the finest mesh that has a uniformly bounded (relative) condition number 

;<(m(^)" a(^)) =0(1), 

and an optimal computational complexity, that is, linear in the number of degrees of freedom A^^ at the 
finest mesh (grid). In order to achieve this goal hieraichical basis methods can be combined with various 
types of stabilization techniques. 

One particular purely algebraic stabilization technique is the so-called Algebraic Multi-Level Iteration 
(AMLI) method, where a specially constructed matrix polynomial p^'^'^ of degree V( can be employed at 
some (or all) levels £. The AMLI algorithm has been originally introduced and studied in a multiplicative 
foi-m, see EH. 

We have the following two-level hierarchical basis representation at level i 



(15) 



i(0 



.^21 



1(0 

^12 



^21 '^22. 

Starting at level (associated with the coarsest mesh), on which a complete LU factorization of the 
matrix A^"^ is performed, we define 

(16) m(0):=a(''). 

Given the preconditioner M^^^^^ at level € — I, the preconditioner M^'^^ at level £ is then defined by 

(17) mW:=lW[/W, 
where 



(18) 



'(0 



At) 

1(0 
/I21 







c 



(0 

22 



u 



(0 



Here is a preconditioner for the pivot block A^^ , and 



(19) 



is an approximation to the Schur complement S = A^'- ^) 



^21 ""ll 



A\^, where A^^-i) 



^22^ is 



the stiffness matrix at the coarse level ( 
satisfying the condition 

(20) 0^p(^)(f)<l, 



1, and p^'-^ is a certain stabilization polynomial of degree Vf 



0<?^1, p(^)(0) = l. 
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It is easily seen that (1% is equivalent to 

(21) C;_J"' =m(^-')"'^W(a(^-i)m(^-i)"\ 



'22 

where the polynomial g(^) is given by 



(22) 



X 

We note that the multilevel preconditioner defined via ([TT] ) is getting close to a two-level method when 
q^^\x) closely approximates 1/x, in which case * A^'-^^^ \ In order to construct an efficient 

multilevel method the action of on an ai^bitrary vector should be much cheaper to compute (in 

terms of the number of arithmetic operations) than the action of A^'^^') \ Optimal order solution algo- 

([) ^ ' 

rithms typically require that the arithmetic work for one application of C22 is of the order 0{Nc-\) 
where Nc-\ denotes the number of unknowns at level £ — I. 

It is well known from the theory introduced in ||3] 5] that a properly shifted and scaled Chebyshev 

polynomial p^'-^ := py^ of degree V( can be used to stabilize the condition number of M^^^ A^'-^ (and 
thus obtain optimal order computational complexity). Other polynomials such as the best polynomial 
approximation of 1/x in uniform norm also qualify for stabilization, see, e.g., Il23l . Alternatively, in 
the nonlinear AMLI method, see, e.g., ||6l, a few inner flexible conjugate gradient (FCG) type iterations 
(for the FCG algorithm see, e.g., Il25l ) are performed in order to improve (or freeze) the residual re- 
duction factor of the outer FCG iteration. In general, the resulting nonlinear (variable step) multilevel 
preconditioning method is equally efficient, and, because its realization does not rely on any spectral 
bounds, is easier to implement than the linear AMLI method (based on a stabilization polynomial). For 
a convergence analysis of nonlinear AMLI see, e.g., |[T8ll26ll27l[3TI . 

Typically, the iterative solution process is of optimal order of computational complexity if the degree 
Vf = V of the matrix polynomial (or alternatively, the number of inner iterations for nonlinear AMLI) at 
level £ satisfies the optimality condition 



(23) 1/V(1-^') <y<r, 

where t % = Nf/Ne-i denotes the reduction factor of the number of degrees of freedom (DOF), and 
y denotes the constant in the strengthened Cauchy-Bunyakowski-Schwarz (CBS) inequality. The value 
of T is approximately 4 in case of the sequence of non-nested N° spaces, which we will construct in the 
next subsection. For a more detailed discussion of AMLI methods, including implementation issues see, 

e.g., mm. 

Remark 3.1. The preconditioner defined in ([TT] ) is of multiplicative form. The introduction of AMLI 
methods was based on the multiplicative form, see E ID |5l , and is commonly used in practice. How- 
ever, it is also possible to choose the preconditioner in the additive form, which is defined as follows 



(24) Aff 



Ci? 

eg 



In this case the optimal order of computational complexity demands that the matrix polynomial degree 
(or the number of inner iterations of nonlinear AMLI) satisfy the following relation 



(25) ^(l+y)/(l-r) < y < T. 

3.2. Hierarchical splitting of the lowest order Nedelec space. The AMLI methods we are considering 
here, for the solution of are based on a proper splitting of the N° subspace of //(Q, curl). The 
particular two-level HE transfoiTnation that induces this splitting was introduced in the context of linear 
nonconfoiTning (Crouzeix-Raviart) elements in |f8l. It was later studied for quadrilateral rotated bilinear 
(Rannacher-Turek) type elements in [12 1. In the following we will briefly describe this transformation 
and study the related splitting of the lowest order Nedelec space. 

Consider two consecutive discretizations Th (coarse level) and T/, (fine level). Figure [T] illustrates a 
macro-element G (at fine level) obtained from a coarse element by one regular mesh-refinement step. We 
see that in this case the corresponding local (and thus the global) finite element spaces 'Vh and "Vh are 
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Figure 1. Macro-element obtained after one regular mesh-refinement step 



not nested. Let (pc = {^i{x,y)}^^^ be the macro-element vector of the nodal basis functions, and be 
the macro-element stiffness matrix con^esponding to G e T = Th. The global stiffness matrix A/j can be 
written as 



A, = 2 Rg^gRg, 



GeT 



where Rq denotes the natural inclusion (canonical injection) of the 12 x 12 matrix Aq for all G in 
T. Using the local numbering of DOF, as shown in Figure [T] (right picture), a macro-element level 
(local) transformation matrix Jq is constructed based on differences and aggregates of each pair of basis 
functions 0, and (pj that correspond to a macro element edge, i.e.. 



(26) 



1 -1 



1 1 



1 -1 



1 1 



1 -1 



1 1 



1 -1 



1 1 



cf. lfT2l . This transformation defines a two-level hierarchical basis ipc locally, namely, (pc = Jg9g- Then 
the hierarchical two-level macro-element matrix is given by 



^gstRg^gRg- 



^G = Jg^gJq' 

and the related global two-level matrix can be obtained via assembling, i.e.. A/, 
Alternatively, one can compute the matrix A/, via the triple matrix product 

(27) Ah = JAhf, 

where the global transformation matrix J is induced by the local transformations, i.e., 

J\g= Jg, VG e T. 

In other words, global and local transformations are compatible in the sense that restricting / to the DOF 
of any macro-element G we obtain Jg- Now, if we number those DOF that correspond to interior nodes 
of the macro elements first, the global two-level stiffness matrix A/, has the 2x2 block structure 



(28) 



Ah 



An An 
All Ml 



where All corresponds to the "interior unknowns". We follow the "first reduce" (FR) approach, see e.g., 
ll8j[T2l, where these interior unknowns are first eliminated exactly. This static condensation step can be 
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written in the form 
(29) 



Ah 



All 
An B 

-1, 



h 




h 



with the Schur complement B = A22 — A2iA^^ A12. Next, the matrix B is partitioned into 2x2 blocks, 
i.e., 



(30) 



B 



Bn 
B21 



B\2 

B22 



where each row (and each column) of Sh corresponds to a difference of two basis functions associated 
with one macro-element edge (face) and each row (and each column) of B22 corresponds to an aggregate 
of such a pair of basis functions. The matrix B22 at level i then defines the coarse-grid matrix A^'-^^^ in 
the AMLI hierarchy, cf. ([T5] ). This algorithm can be applied recursively on each level ^ = L, L— 1,...,1. 
The resulting algorithm is then of optimal computational complexity, see e.g., |[22l Remark 3.1]. 

3.3. Local analysis. In the two-level framework we denote by 'Vi and 'V2 the subspaces of the finite 
element space 'Vi,. The space 'V2 is spanned by the coarse-space basis functions (aggregates) and "Vi is 
the complement of 'V2 in "V/,, i.e., "V/j is a direct sum of "Vi and 'V2: 

(31) T/, = %®^2. 

A measure for the quality of this splitting is the constant y in the strengthened CBS inequality, which is 
defined by the relation 

'H{u,v) 



r = cos(%,T2) := 



sup 



It is well known (see, e.g., 13) that 7 can be estimated locally over each macro element G, and that 
y = maxG yc^ where 



7g 



sup 



The spaces 'Vi (G), 'V2{G), and the bilinear form Hciu, v) correspond to the restriction of "Vi, 'V2, and 
9{{u, v), respectively, to the macro element G. 

We perform this local analysis on the matrix level, where the splitting (|3TI ) is obtained via the two-level 
hierarchical basis transformation described in Section l3T2l and the space "Vh corresponds to the choice of 
lowest order Nedelec elements. In this setting the upper left block of A/, is block-diagonal (with diagonal 
blocks of All of size 4x4, which can be associated with the interior nodes {1, 2, 3, 4} in the right picture 
of Figure [T|). Therefore, we first compute the local Schur complements arising from static condensation 
of the interior DOF and obtain hereby the (8 x 8) matrices Bq- Next we split each matrix Bq as 

differences 
aggregates ' 

written again in two-by-two block form (with blocks of size 4 x 4). We have thus reduced the problem of 
estimating the CBS constant of the splitting (|3TI ) to a small-sized local problem that involves the matrix 
Bg- Following the general theory, see ||2j[Tl], to estimate the CBS constant y, it suffices to compute the 
minimal eigenvalue of the generaUzed eigenproblem 

(32) 'S'gVg = /lG,min5G,22VG, VVg, 

where So = ^0,22 — Bq 2iB~^^^Bg,i2- The CBS constant y can then be estimated as follows 

(33) y^ ^ maxy^ = max(l - /icmin)- 



Br 



Bgm 

Bg,2\ 



Bgm 

Bg,22 



max(l 

GeT 



Note that the matrix Bg,ii is a well conditioned matrix, see Figure |21 and therefore, it can be inverted 
cheaply, either by an iterative process or by, for example, an incomplete LU factorization OUII . which is 
denoted by B'^^ in Figure |21 

We now first prove two auxiliary (stand-alone) results on algebraic sequences, which we will use to 
bound the CBS constant y. 
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_K(B„) 

_K((B' r'B,. 



2.5- 



6 

Number of levels 



Figure 2. Condition number of the matrix Bq h 



Lemma 3.2. For all e > 0, consider the coupled sequences 

(34a) bQ = e -6, ao = 2e + 6 = 2{bo + 9), 

(34b) be+i = -b]/ai, a^+i = 2ac + be+u € = 0,1,2,.... 

Let r[ = b[/ai. Then we have 

(35) bc+i/ae = -r], ae+i/ac = 2-r], rt+i = -r]/{2 - r]). 

Moreover, the following bound holds for all £ = 0,\,2, . . . 

(36a) ai > . . .a\ > a^ > 6, 

(36b) ^ ^ ... ^ ^ < 1. 

Proof. Using the definition of r( in (I34bl ). we get b(+\/a( = —r^, and thus a{+\/a[ = 2 — r^. The last 
relation of (l35T l then immediately follows. From (I34bl i and (1351 ) we also note the following relations 

(37) af+i - be+i = 2ae, ai+y + bf+y = —{a] - b]) = 2ae{l - r^,,), = i _ ^2 

Clearly, for e > 0, we have a^ > 6, and since ro = bo/ao = {e — €)/{2e + 6), it is easy to see that 
— 1 < ro < 1/2. The latter also implies that ^ < 1. We now prove the remaining bounds using 
induction. 

^ = 0. Since a\/aQ = 2 — > 1, we have ai > > 6. Moreover, r\ = —r^/{2 — r^). This implies 
that — 1 < ri ^ 0, and thus ^ Tj < 1. Furthemiore, when ro t 0, we have 

And, since r\ = if ro = 0, we have r^ ^ r^ < 1 . 
i = n. Assume that the relations (|36] | hold for t = n. Since an+\/an = 2 — r^ > 1, we have a,j+i > > 
6. Moreover, r„+i = — r,^/(2 — r,^). This implies that —1 < r^+i ^ 0, and thus ^ r^^^ < 1. 
Also, when r„ i= 0, we have 

2 ^ I 'n \ n+l ^ 'n , 

'"^'-\2-rl) rl ~ {2-rlY ' 

And, since r„+i = if r„ = 0, we have r^^^j ^ r^ < 1. 
This concludes the proof. □ 

From Lemma 1221 we also note the following bounds 

(38) - 1 < ro < 1/2, and - 1 < rf ^ V ^ = 1, 2, . . . . 
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Lemma 3.3. Let e > and the sequences ac and be be as defined in Lemma \J?2\ Then for 

(39) c] = ^^(^^ + ^^^ 

{a] - 36) (a^ - Zj^)' 

the following bounds hold for all £ = 0, 1,2, .. . 

(40) < c^_j < ...<c]<cl< 3/8. 
Proof From ao = 2e + 6 and bo = e — 6, we have ao — bo = e + 12, ao + bo = 3e, ao — 6 = 2e, and 

.2 
-0' 

27 

(f' + 6)(7+ 12) 
Now 

36 ((fli + Z7i)(fl'f) - 36)(ao - ^o) - («o + ^o)(ai - 36)(ai - Z^i)) 

^ — ^ — 

(a^ - 36) (ai - bi){al - 36) (ao - ^o) 
Substituting the values of ao, ai, ^o and bi, and after some lengthy, but simple calculations, we find that 



ao + 6 = 2{e + 6). Substituting these relations in the definition of c^, we get 
(41) (^l = , Jl , < 3/8. 



^2 2 



108f' (-9^2(312 + 80^ + 5^2)) 
e + 3)(aj - 36) (ai — bi){al — 36) (ao — bo) 

2 „2 



Since the denominator is a positive quantity, we get — Cq < 0, and thus 
(42) < 3/8. 

For remaining bounds, we again use induction. Note that, using (l37l) we get 



(43) 



36{ae+i + be+i) 36(1 



r 



.2^ 



(a^^j - 36)(a^+i - b(+i) {aj_^^ - 36) 
Therefore, to show that c^^j < 3/8, it suffices to show that 

(44) a^^i -36 > 96(1 -r^). 

Since Cj < 3/8, we clearly have aj — 36 > 96(1 — r^). Now assume that the relation (|44] | holds for 
{ = n- I, i.e., 

(45) al - 36 > 96(1 - J. 
Multiplying (|45] ) by (2 — ^2)2 and subtracting 36 from both sides we get 

(2 - r2)2a2 _ 36 > 36(2 _ ,2)2 ^ 96(1 _ ^2_^)(2 - r^)^ - 36 

al^, - 36 > 96 ((2 - rl)\n/S - r,ti) - 3/8) . 
We need to show that (2 - rfy{ll/S - rl^) - 3/8 > 1 - rl, i.e., 

(46) gn := (2 - r2)2(ll/8 - r^j) + - 11/8 > 0. 
From the recurrence relation on r„ from (l35T l. we have 

Substituting these relations in g„, and after some lengthy calculations we obtain 

Now for j e [0, 1), we have 

l-rti>0, 2-rti>0, 66 - 64rt, > 0, 15r^i - r,ti ^ 0, 

which proves that g„ > 0, and that a^^^ — 36 > 96(1 — r^). Therefore, the inequality (|44] | holds for all 
^ = 0,1,.... 
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To prove the mono tonicity of c^, we show that 
(48) fc := c]^Jc] < 1. 

Using (l43l) we get 

{\-r]){a]-3€) 



fe 



;i-^ti)K+i-36) 
Multiplying numerator and denominator by (2 — r^)^, we obtain 

(1 - r^) ((2 - r]fa] - 36(2 - r]Y) 



fi 



(1_^) (a^+,- 36 + 36(l-(2-r^)^) 

(l-'-Li) 36)(2-r2)2 

(l-'^) 36(l-r2)(l-(2-r2)2) 



(1 - J(2 - r]Y (1 - r2_ J(a2^, - 36)(2 - r]Y ' 



Now since c]^^ < 3/8, we have (1 - rj)/{aj_^_^ - 36) < 1/96 from dSll. Therefore, 



[I - r2) 36(1 - (2 - r2^2^ 

2 \(n _ ^2\2 QAM _ »-2 VO _ ,-2^2 



;i-r2_J(2-r2)2 96(l-r2_j)(2 

+ (2-^2)2) ll/8-r2-^(2-r?)2 



(l-r2_J(2-r2)2 (i_,2_^)(2-,2)2 



This gives 



ll/8-r2-^(2-r2)2-(l-.2_J(2_,2)2 



Using (l46l) we therefore get 



ll/8-r2 + (2-r2)2(-ll/8+r2_^) 
(l-r2_^)(2-r2)2 



since > 0, 1 — j > 0, and (2 — r^)^ > 0. This proves (|48] ) and concludes the proof. □ 

The sequences , , rc and are plotted in Figure |3] We aie now in a position to prove the following 
theorem which provides a theoretical estimate that holds on all levels of recursive splitting of the 
subspace of H{D., curl). 

Theorem 3.4. Consider the bilinear form ([T|) where < a,p < co, and the related discrete prob- 
lem (O on the N" subspace of H{Q., cm\) where the underlying rectangular mesh is uniform. Under 
these assumptions the CBS constant y related to the hierarchical splitting diil ) has the upper bound 

7 ^ 7g < -^3/8, which holds for each step of the recursive hierarchical splitting. Moreover, 

monotonically strictly decreasing and has an upper bound of \f^ij%for all i = 0,1,. .. ,L, i.e., 

(49) y«') < r(') < . . . < r^') < . . . < 7^"^-'^ < y(^^ < 

Proof. In order to prove this uniform bound for y we study the generalized eigenproblem (|32] ). At level L 
of the finest discretization the macro-element matrix Ac, which is the same for all G in Th^ for a uniform 
mesh, can be represented in the form 

(50) 4''^=./g( 2 RIa^^'^Rk]^ 
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Figure 3. a^, b(, rc, and for e = lO™", where niQ = {2, 1, 0, ... — 11, —12} (left to right) 



where 



(51) 



ao bo —6 6 

^(L) ^ _^ bo ao 6 -6 
^ 6/i2 -6 6 £?() ^0 

6 —6 bo ao 

with ao = 2e + 6, Zjq = e — 6, e = a-Zi^ > 0, with k = a/p, and the local transformation matrix Jq is 
defined according to (l26l ). 

The lower-right 4x4 block of the matrix Bq and the Schur complement Sq for the first splitting (at 
level L) are to be found 
(52) 



B 



G,22 



6h^ 



PO 


qo 


-3/2 


3/2 




PO 


3/2 


-3/2 


3/2 


3/2 


Po 


qo 


3/2 


-3/2 


qo 


Po 



,(L) 



6/z2 



■^0 

to 



to 



-3/2 



-3/2 3/2 
3/2 -3/2 



^0 3/2 
So 



to 



3/2 
-3/2 
to 
So 



with 



Po = ao/2 - bl/4ao, 

108c?o - 2al + llbo + aobl 



So 



qo = -bl/Aao, 

36(30 + 12bo + aob^ 



144 _ 



to 



144 _ 



The generalized eigenproblem (l32l ) has two different two-fold eigenvalues, namely Ai2 = l and 

aoial - aobo - 72) 



which shows that 
(53) 



(aq - 36) (ao - ^o) ' 

36(ao + ^o) 



{al - 36) (ao - ^o) 



Note that the coefficient yS does not appear in the bound for y since the factor ^ appear in both the 
matrices of the generahzed eigenproblem (|32] |. and thus does not affect the eigenvalues. 
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Now in Older to compute a similar bound for the second splitting (at level L — 1) we have to use the 
relation := B^22- general, for the {£ + 1)* splitting (at level L — () the relation 

(54) 4-'^ := 

is to be used in the assembly of A^^, i.e., 



(55) 4"'-.'= s 4^ 



Repeating the computations, we find that the relation (155] ) holds for all levels ^ = 1, 2, . . . , L — 1, L, and 



the element stiffness matrix A„ (after £ coarsening steps) is given by 



(56) 4-'^ - 



ac b( —b b 

be a[ 6 —6 

—6 6 a[ b{ 

6 —6 b( ac 



e G, VG c Th, 



where the sequences af and b( are defined in (l34l ). Thus, the bound for yc at level L — £ reads 



^^"^^ ^ («2-36)(a, 

The result (|49] l then follows by taking 7^ = in Lemma [ 

Remark 3.5. The curves in the lower-right part of Figure |3]show the behavior of 7^ (defined by (|57])). 
We observe that 7^ approaches zero when the splitting is applied many times (increasing £ from left to 
right), which means that the two subspaces "Vi and 'V2 in (|3TI ) become increasingly orthogonal to each 
other as the recursion proceeds. Therefore, on (very) coarse levels, the upper bound 3/8 for 7^, and thus 
for 7, is quite pessimistic. 

4. Numerical results 

All the numerical experiments presented in this section are performed using Matlab R2012a on an HP 
Z400 workstation with 6 core 3.33GHz CPU and 12 GB RAM. For all the numerical experiments, we 
consider a mesh of square elements of size h = 1/8, 1/64, . . . , 1/2048 (i.e., up to 8, 392, 704 DOF for 
the finest level). We use a direct solver on the coarsest grid that consists of 2x2 elements. Hence, the 
multilevel procedure is based on 1 to 9 levels of regular mesh refinement (resulting in a ^-level method, 

= 3, . . . , 1 1). The initial guess is chosen as a zero vector, and the stopping criteria is chosen as 



(||r("-')||/||r(°)||) ^ 10- 



where nn is the number of iterations reported in the tables. The average residual reduction factor p, which 
is also reported in the tables, is defined as 

p:= f||r(«-')||/||rW| 



Example 4.1. Consider the model problem ([T]) in a unit square, and fix the coefficients a = fS = 1. The 
problem data is chosen such that the exact solution is given by u = (;r sin /ta; cos ;ry, — ;7rcos;r.xsin;rj)^. 

For the W-cycle method, we chose two-types of stabilization polynomials q^'^K One is based on 
Chebyshev polynomials (see, e.g., ||4j[19j|3Tl, denoted in the tables by T), for which the polynomial q^^^ 
is defined as 2/ {s — b) — t/{s — b)^, where s = -y/l + b + b^ — y^^, and b is some constant estimating 
the upper bound of the condition number of preconditioned Bn block. The other one is based on the 
polynomial of best uniform approximation to l/;c (see, e.g., 1231 . denoted in the tables by X), for which 
the polynomial q^^'' is defined as (2 — 7^)/(l — 7^) — t/{l — 7^). The results for the V-cycle and W- 
cycle multiplicative AMLI method are presented in Table [T] The second column confirms the error 
convergence behavior. We see that for decreasing h the growth in the iteration number for V-cycle is 
moderate (as expected), whereas both the W-cycle versions (T and X) exhibit /2-independence. Moreover, 
the total time (factorization and solver) reported in eighth and eleventh columns also confirms that both 
the versions of W-cycle are of practical optimal complexity (slight increase in time may be attributed to 



ALGEBRAIC MULTILEVEL PRECONDITIONING IN H(n. curl ) 13 

the implementation issues). We note that the X version W-cycle gives slightly better results than the T 
version W-cycle. 



Table 1 . Convergence results for multiplicative AMLI, a = /3 = I, x = u — uj, 







V-cycle 


W-cycle (T) 


W-cycle (X) 


1/h 


curl y r^/'n^ 


"it 


P 


^sec 


"it 


P 


^sec 


nit 


P 


^sec 


8 


0.15946423 


7 


0.049 


0.00 


7 


0.049 


0.00 


1 


0.049 


0.00 


16 


0.08005229 


8 


0.094 


0.01 


8 


0.083 


0.01 


8 


0.094 


0.01 


32 


0.04006629 


10 


0.143 


0.01 


9 


0.104 


0.01 


8 


0.097 


0.01 


64 


0.02003817 


11 


0.174 


0.04 


9 


0.105 


0.05 


8 


0.100 


0.04 


128 


0.01001971 


12 


0.201 


0.14 


9 


0.108 


0.16 


8 


0.095 


0.14 


256 


0.00500993 


13 


0.224 


0.54 


9 


0.109 


0.55 


8 


0.088 


0.51 


512 


0.00250498 


14 


0.246 


2.41 


9 


0.110 


2.22 


8 


0.083 


2.09 


1024 


0.00125249 


14 


0.267 


10.74 


9 


0.110 


9.35 


8 


0.078 


8.99 


2048 


0.00062624 


16 


0.313 


49.79 


9 


0.110 


40.29 


8 


0.073 


38.74 



We now test the AMLI method with additive preconditioning. The results for the V-cycle and both 
the W-cycle additive AMLI methods are presented in Table |2] We also present the results for nonlinear 
variant of AMLI method, see, e.g., 121 IH [HI [19l |25l HH |27l, in the last three columns (denoted in the 
tables by A'^, W-cycle referring to two inner iterations). Surprisingly, in the additive form, the T version 
W-cycle gives much better results than the X version W-cycle, where the latter appears to be stabilizing 
only towards the very fine mesh (many recursive levels). This can be attributed to the fact that for the 
additive preconditioning, for the choice of y = -^3/8, we require that v > ^ {\ + 7)/(l — j) > 2. The 
results of nonlinear W-cycle further improve the results of T version W-cycle (linear). Since the nonlinear 
W-cycle AMLI method gives the best results (and is free from parameters b and y), in the remaining 
numerical experiments we will only present the results from nonlinear W-cycle AMLI method. 



Table 2. Convergence results for additive AMLI, a = = \ 





V-cycle 


W-cycle {T) 


W-cycle (X) 


W-cycle (AO 


l/h 


"it 


P 


^sec 


"it 


P 


^sec 


"it 


P 


^sec 


"it 


P 


^sec 


8 


10 


0.153 


0.00 


10 


0.153 


0.00 


10 


0.153 


0.00 


10 


0.153 


0.00 


16 


17 


0.300 


0.01 


17 


0.299 


0.01 


17 


0.299 


0.01 


12 


0.208 


0.01 


32 


20 


0.391 


0.02 


19 


0.346 


0.03 


23 


0.446 


0.03 


12 


0.209 


0.03 


64 


25 


0.472 


0.06 


19 


0.372 


0.08 


31 


0.550 


0.13 


12 


0.197 


0.08 


128 


30 


0.538 


0.21 


21 


0.386 


0.26 


44 


0.653 


0.47 


11 


0.179 


0.23 


256 


34 


0.575 


0.87 


19 


0.377 


0.79 


56 


0.712 


1.82 


11 


0.167 


0.76 


512 


39 


0.617 


3.99 


19 


0.361 


3.04 


60 


0.735 


6.85 


9 


0.127 


2.61 


1024 


44 


0.657 


19.20 


19 


0.362 


12.46 


65 


0.751 


28.90 


9 


0.117 


10.65 


2048 


50 


0.685 


91.72 


19 


0.371 


52.99 


65 


0.752 


121.49 


8 


0.098 


43.03 



Example 4.2. Consider the model problem ([T|l in a unit square, fix the coefficient jS = \ and take 
a = 10"^ for mo = {—6, —3, 0, 3, 6}. The right hand side (RHS) vector is all ones. 

The results for the multiplicative AMLI method for varying a are presented in Table |3] for V- and 
nonlinear W-cycle. We see that the V-cycle shows some effect of a, with a moderate growth in the 
number of iterations for decreasing h, however, the nonlinear W-cycle is independent of h, and is fully 
robust with respect to a. Note that towards very large values of a, the system matrix is well-conditioned, 
and the hierarchical splitting approaches orthogonal decomposition, therefore, the V-cycle method also 
exhibits optimal order complexity. 

Example 4.3. Consider the model problem ([T|) in a unit square, fix the coefficient a = \ and take 
P = lO'"o for mo = {-6,-3,0,3,6}. The RHS vector is all ones. 
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Table 3. Convergence results for multiplicative AMLI, jS = 1 





"it 


a — > 


10 


-b 


10 


-i 


10 


I) 


10 


i 


lO'' 


\/h 


V 


w 


V 


W 


V 


W 


V 


W 


V w 


8 


9 


9 


9 


9 


9 


9 


4 


4 


2 2 


16 


12 


10 


12 


10 


12 


10 


7 


6 


2 2 


32 


15 


10 


15 


10 


14 


10 


9 


8 


2 2 


64 


17 


10 


17 


10 


16 


10 


11 


9 


2 2 


128 


20 


9 


20 


9 


17 


9 


12 


9 


3 3 


256 


22 


9 


22 


9 


18 


9 


14 


9 


4 4 


512 


26 


9 


26 


9 


21 


9 


16 


9 


6 6 


1024 


28 


9 


28 


9 


23 


9 


17 


9 


8 8 


2048 


28 


9 


31 


8 


25 


8 


20 


8 


10 8 



Table 4. Convergence results for multiplicative AMLI, a = 1 





"it 




10 


-b 


10 


-i 


10 


u 


10 


i 


10" 


1//2 


V 


w 


V 


W 


V 


W 


V 


W 


V w 


8 


2 


2 


4 


4 


9 


9 


9 


9 


9 9 


16 


2 


2 


7 


6 


12 


10 


12 


10 


12 10 


32 


2 


2 


9 


8 


14 


10 


15 


10 


15 10 


64 


2 


2 


11 


9 


16 


10 


17 


10 


17 10 


128 


3 


3 


12 


9 


17 


9 


20 


9 


20 9 


256 


4 


4 


14 


9 


18 


9 


22 


9 


22 9 


512 


6 


6 


16 


9 


21 


9 


26 


9 


26 9 


1024 


8 


8 


17 


9 


23 


9 


28 


9 


28 9 


2048 


10 


8 


20 


8 


25 


8 


31 


8 


28 9 



The results for the multiplicative AMLI method for vaiying y6 are presented in Table |4] for V- and 
W-cycles. The results are qualitatively the same as in Table [3] for varying a, with the parameter value 
reversing the behavior of the solver. 

Example 4.4. Consider the model problem ([1]) in a unit square, and fix the coefficient /? = 1. The coeffi- 
cient a is chosen as 1 in [0, 0.5]^ |J(0.5, 1]^ and k elsewhere, where k = 10'"", andmo = {—6, —4, —2, 0}. 
The RHS vector is all ones. 

Finally, the results for the multiplicative AMLI method for the case with jump in the coefficients, 
which are presented in Table |5]for V- and nonlinear W-cycles, show robustness with respect to the jump 
in the coefficients. 

5. Conclusion 

We have presented an optimal order AMLI method for problems in the H{Q., curl ) space. The main 
result of our local analysis (Theorem I3.4l i shows that a second order stabilization polynomial (or two 
inner iterations in nonlinear method), i.e., a W— cycle, is sufficient to stabilize the AMLI process. We 
derived explicit bounds for the multilevel behavior of y that are robust with respect to the coefficients 
in the H{Q., curl) model problem. The presented numerical results (including jumping coefficients case) 
confirm the robustness and efficiency of the proposed method. 

Acknowledgements. The author is very grateful to Dr. Johannes Kraus (RICAM, Linz) for insightful 
discussions on AMLI methods. 
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